Kinetic theory of quantum transport at the nanoscale 
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We present a quantum-kinetic scheme for the calculation of non-equilibrium transport properties 
in nanoscale systems. The approach is based on a Liouville-master equation for a reduced density 
operator and represents a generalization of the well-known Boltzmann kinetic equation. The system, 
subject to an external electromotive force, is described using periodic boundary conditions. We 
demonstrate the feasibility of the approach by applying it to a double-barrier resonant tunneling 
structure. 
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I. INTRODUCTION 

There is currently large interest in the development 
of electronic devices that operate at the molecular level, 
since these may revolutionize the electronic industry 1 . 
Modeling electron transport at the nanoscale poses great 
theoretical challenges because in this regime the semi- 
classical Boltzmann kinetic equation commonly used to 
describe transport is not applicable. Boltzmann's equa- 
tion treats the electrons in terms of their classical prob- 
ability distribution in phase space, and deals with colli- 
sion events quantum-mechanically using Fermi's golden 
rule^. This formalism explains why conduction electrons 
subject to a uniform electric field do not accelerate in- 
definitely, but settle in a steady state regime of constant 
current, as a consequence of inelastic collisions with the 
lattice that lead to energy dissipation. 

To describe transport in devices having spatial di- 
mensions of the order of the electron wavelength, Boltz- 
mann's equation should be replaced by a quantum me- 
chanical Liouville-master equation for a reduced density 
operator^^^. Such quantum kinetic approaches to 
electron transport have to deal with several difficulties. 
Probably most severe are those related to the description 
of boundary conditions between the device active region 
and the contacts: In most common approaches electrons 
can be exchanged between the device and its environment 
through left (L) and right (R) contacts with the elec- 
trodes of a battery. The L(R) injected electrons are as- 
sumed to have thermal equilibrium distributions defined 
by the temperature T and the electrochemical potentials 
/ i L(fl) °f the electrodes. In self-consistent calculations, 
however, modifications of the distribution of the injected 
carriers, like drifting the momentum of the distribution 
or displacing the chemical potentials of the contacts, are 
necessary to avoid unphysical charging of the contacts^. 
Overall the issue of contacts and boundary conditions is 
crucial in all models of quantum transporliSiZ&SiiSiiiii^. 

In this paper we present a kinetic formulation which 



avoids these difficulties by adopting periodic boundary 
conditions. The active region together with a part of the 
contacts is placed in a unit cell that is repeated period- 
ically. One can imagine such a geometry as an infinite 
ring, without open boundaries. In this way we deal with 
a closed circuit with no exchanges of electrons with the 
environmentpi&iiiiii. A current is induced in the circuit 
by an externally applied electromotive force. 

In our approach coupling between the electrons and 
the lattice is included explicitly to prevent the electrons 
from accelerating indefinitely. This coupling is treated in 
perturbation theory (weak coupling). Furthermore, the 
time scale of inelastic electron-phonon collisions is typ- 
ically much faster than the current relaxation time and 
is eliminated by a coarse graining procedure. The re- 
sulting Markovian dynamics of the reduced electron den- 
sity operator is described by a Liouville-master equation, 
which effectively generalizes Boltzmann's equation at the 
nanoscale^. We treat the diagonal and off-diagonal ele- 
ments of the density operator on equal footing. 

One common problem of kinetic master equations that 
describe collisions using Fermi's golden rule is that cur- 
rent continuity is violated in inhomogeneous media when 
the current density is defined in the standard way*S. In 
our approach, continuity is exactly satisfied, by prop- 
erly taking into account the current contribution due to 
electron-phonon collisions, as described elsewhere^. 

The reminder of this paper is organized as follows. 
In section [H] we describe the methods that we use. In 
particular, we recall how a Liouville master equation for 
an electron in contact with an external heat bath is de- 
rived (section III Afl . We then discuss electric fields (sec- 
tion ^TOJ, the definition of currents (section III C|l . and 
a generalization of the formalism to many particle sys- 
tems Ijll D0 . Finally, in section ITTT1 we demonstrate the 
scheme with a numerical application to a system con- 
sisting of a one-dimensional wire connected to a double- 
barrier resonant tunneling structure (PBRTS). This ex- 
ample illustrates the key features of our approach and 
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the differences between our treatment and open bound- 
ary formulations£i2ii2iii. 



II. METHOD 

A. Liouville master equation 

For simplicity we use in the following a one- 
dimensional notation, but extension to three-dimensions 
is straightforward. A ring of length L is a simple peri- 
odic structure. Ideally we should consider an infinite ring 
(L —y oo). In numerical applications, however, L must be 
finite. This could be achieved by placing many periodi- 
cally repeated unit cells of length L in a ring, but in the 
following we shall consider only the case of a single peri- 
odic unit having length equal to the perimeter of the ring. 
In absence of the external electric field the single-electron 
Hamiltonian is H = l/2p 2 + U(x) = X^n e " c n c i' where 
U(x) is a potential with periodicity L, the sum is over 
eigenstates of energy e„, and c4(c„) are electron creation 
(annihilation) operators. Atomic units (e = h = m = 1) 
are adopted throughout. The harmonic phonon bath has 
Hamiltonian R = ^ Q o; Q aJ 1 ,a ct , where the sum is over 
eigenmodes of frequency uj a , and a],(a a ) are phonon 
creation (annihilation) operators. The electron-phonon 
coupling potential is V = Ero.n.a 7m,n c n c m a a + h - c -i in 
terms of matrix elements 7^ n {h.c. denotes the hermi- 
tian conjugate). The Hamiltonian of the total system 
(electron + bath) is Ht = Hq + R + V. The density op- 
erator for the total system £ satisfies the Liouville equa- 
tion — — i [Ht,^]- The reduced density operator for 
the electronic subsystem, S, is obtained from £ by trac- 
ing out the bath degrees of freedom, i.e. S — Tr^E. 
The equation of motion for S can be derived using a 
standard procedur e) 18 ! 19 in which V is treated to second 
order in perturbation theory, and a Markov approxima- 
tion is made, whereby electron-phonon collisions become 
instantaneous processes described by Fermi's golden rule. 
The bath is assumed to be in thermal equilibrium at all 
times, i.e. any retroaction of the electron system on the 
bath is neglected. Then the Liouville-master equation for 
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^- = -i[H ,S]+C[S] 
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(1) 



Here the term that describes collisions with the bath 



is 



m,n 
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where L nm = c ra c„ 
are given by: 



and the transition probabilities T r 
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Here n w is the thermal occupation and T>(uS) is the 
density of states of phonons with frequency to. Eq. © 
ensures detailed balance leading the electron subsystem 
to thermal equilibrium in absence of an applied external 
electric field. 

In our scheme, we calculate the time evolution of the 
density matrix, including diagonal and off-diagonal ele- 
ments. In approaches using open boundary conditions, it 
is common practice to limit the description to the diag- 
onal elements of the density matrix onlyi. It should be 
stressed, however, that the representation of the density 
matrix is very different when open boundary conditions 
are used instead of periodic boundary conditions. In open 
boundary approaches the different chemical potentials of 
the L and R contacts cause a beaking of symmetry be- 
tween the contacts, and the final state is largely deter- 
mined by the different occupations of the L and R inci- 
dent channels. These occupations are given by the diago- 
nal elements of the density matrix in the open boundary 
representation. In our approach, where periodic bound- 
ary conditions are adopted, we need to use unperturbed 
eigen] 'unctions as basis states. These states do not carry 
current, and the final, current-carrying state of the sys- 
tem is dominated by the off-diagonal elements of the den- 
sity matrix in this basis. The computational complexity 
remains limited, however, because not all the elements of 
the density matrix need to be considered, but only those 
that fall within a finite energy window around the Fermi 
energy. 



B. External static electric field 

We include an external uniform static electric field £ 
by letting p — > p — l/cA(t). The Hamiltonian becomes 
H(t) = 1/2 (p- l/cA(t)f+U{x). The vector potential is 
A(t) = —c£t. This choice of the gauge (v-gauge) is com- 
patible with periodic boundary conditions. We avoid the 
difficulty with a Hamiltonian that grows indefinitely with 
time by systematically performing gauge transformations 
in which the electronic wavefunctions are multiplied by 
a phase factor exp(— iA(t)/c ■ x) and A(t) in the Hamil- 
tonian is restored to its initial value A(t = 0)2£. In this 
way all the time dependence is transferred to S and the 
Hamiltonian remains constant. In order not to violate 
periodic boundary conditions the gauge transformations 
are only possible at times that are integer multiple of 
T£ = jg. The time T£ is arbitrarily small for L — > oo, 
but is always finite in numerical applications. In practice 
we proceed as follows. Let A(0) = at t = 0. First we 
integrate Eq. for the time period T£ using the Hamil- 
tonian propagator with the electric field, which corre- 
sponds to the first term in the right hand side (r.h.s.) of 
Eq. (JIJ. Then we apply a gauge transformation to set 
the vector potential to zero. Finally, we propagate S for 
the same period T£ using the collision propagator, which 
corresponds to the second term on the r.h.s. of Eq.JIJ. 
The result is a density matrix S(ts). This procedure is 
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repeated at all subsequent time steps in increments of 
T£. In the limit of Tg — > this procedure converges to 
the exact kinetic propagation. 

Our approach in which a static external electric field 
acts as driving force for the electrons, is well suited to 
dealing with systems where the voltage drop occurs along 
one given direction in space. Complex geometries in de- 
vices with many-terminals 8 may not be tractable in this 
framework. However, many current applications of nano- 
scale transport, like e.g. break junctions, STM tips, or in- 
terlayer tunneling, can be represented within our scheme. 



C. Hamiltonian and dissipative currents 



The independent particle system is described by a 
single particle density operator S , given by S 1 = 
J2im flm I0( m l m tne basis of the eigenstates of H , i.e. 
the equilibrium single particle eigenstates. The expan- 
sion coefficients //,„ are related to the many particle oper- 
ator S via firp, = Tr [Sc^c/] . This procedure is discussed 
e.g. in Ref. l2ll The equation of motion for // m then 
follows from Eq. (HJ). Using Tr [Sc^c^CpCq] = f qn f P m 
(neglecting the exchange term —f qm fpn), one finds: 



dfn 



dt 



(Sr, 



fnm) ^ (-Tnp "1" fmp) fpp (6) 



The current density associated to Hamiltonian propa- 
gation is 



j(x;t) = Tr [S(t)J(x;t) 
1 



J(x;t) = -[(p-l/cA(t))8(x-x) 
+S(x- x)(p-l/cA(t)) 



(4) 



where the current (as well as all other observables) is 
averaged over the period t$ . It is important to note that 
these definitions do not yield a current that satisfies the 
continuity equation 



dS(t) 



dt 



(5) 



In contrast, Eq. © is satisfied by the physical current 
jr{x; t) = j(x; t) +jc(x; t), which includes a contribution 
jc due to the coupling with the phonon bath. An explicit 
expression for jc{ x \ t) is given in Ref.ll6l In all numerical 
applications below we report the total physical current 
3T(x]t). 



D. Many particle formulation 

In the many electron case, Eq. is still valid, but 
should be reduced to an effective single-particle form 
in practical applications. This is straightforward if we 
adopt the time dependent Hartree approximation. In 
this approximation the many-particle system is mapped 
onto a system of non-interacting particles, subject to a 
time-dependent Hartree potential defined in terms of the 
instantaneous electronic charge density^. The Hartree 
interaction is crucial to modeling charge rearrangement 
effects when the system is driven away from equilibrium. 

The single-particle Hamiltonian in presence of an elec- 
tric field is H(t) = 1/2 (p- l/cA(t)) 2 + U(x) + A$(x;t), 
where A$ is the periodic Hartree potential given by 



A$(x;i) = / dx 



I An(x' \t) 



where An is the charge den- 



fnm (r 
p 



pn 



+ Tpm) (1 fpp) 



sity perturbation induced by the field. 



Here H nm (t) are Hamiltonian matrix elements between 
the equilibrium single-particle eigenstates. In Eq. © the 
diagonal elements f nn satisfy < f nn < 1 in accordance 
with Pauli's principle, and the total number of electrons 
is conserved. At equilibrium and in absence of an electric 
field, Eq. © gives f„ m = S nm / (1 + exp (/3(e„ - /i))), i.e. 
a Fermi-Dirac distribution in which the chemical poten- 
tial n is determined by the condition J2 n fnn = ^Vj where 
N is the total number of electrons. In presence of a static 
electric field Eq. @) leads to a stationary distribution in 
which the effect of the field is balanced by the effect of 
inelastic phonon scattering. In the following wc find the 
stationary distribution by setting / = in Eq. JHJ and 
solve the resulting equation by a self-consistent proce- 
dure. 



III. APPLICATION TO A DBRTS 

To illustrate how this works in practice we consider a 
DBRTS as a test system. The potential U(x) in absence 
of an applied electric field is shown in the upper panel of 
Fig. ^ Outside the barrier region, we assume a carrier 
density of 4.3 • 10 18 cm~ 3 , and choose an effective mass of 
0.1m e and a dielectric constant of 10. Due to the large 
dimensions of the DBRTS in the directions perpendic- 
ular to the barriers, an effective one-dimensional model 
can be assumed. We adopt a model in which the phonon 
density of states 'D(uj) oc lu 2 , and the electron-phonon 
couplings "fnm = 7o for all states n, m. We choose values 
of 0.136meV and of 0.2l8meV, respectively, for 70, and a 
temperature of T — 25.3/^, such that ksT is comparable 
to the electronic level spacing at the Fermi energy in our 
finite system (L = 244 nm) . We solve Eq. © for steady 
state behavior, i.e. — 0, and determine the Hartree 

potential A<&(x) self-consistently, for different values of 
the applied external field 8. The total potential acting 
on the electrons when £ corresponds to peak current is 
shown in the middle panel of Fig. The total potential 
includes U(x), the externally applied bias (represented 
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FIG. 1: Upper panel: the DBRTS in the absence of an electric 
field. Middle panel: total potential in the presence of exter- 
nal electric field for low 70 (see text). Lower panel: spatial 
distribution of the current. Short dashed line: Hamiltonian 
current j(x). Long dashed line: Current due to inelastic col- 
lisions jc{x). Solid line: Total current jr(x) = j(x) + jc{x). 
Within the numerical precision of the calculation jr is con- 
stant as required by continuity All plots use a non-linear 
scale for the position to magnify the barrier region. 



by a linearly varying potential 4>e(x) = —£ ■ x), and the 
Hartree potential. Notice that we adopt here, for visu- 
alization purposes, the position gauge (x-gauge) which 
is obviously incompatible with periodic boundary con- 
ditions. As we can see in the figure, the total potential 
has the same qualitative behavior found in self-consistent 
calculations using open boundary conditions for similar 
model systems. In our approach, however, the voltage 
drop across the barrier is an output rather than an input 
of the calculation. Another distinctive feature is that 
we observe not only a potential drop at the tunneling 
structure but also a small voltage drop along the wire. 
The latter is very small (about 5 % of the total in the 
case shown in Fig. ^) , indicating that the external field is 
nearly completely screened inside the wire, as one should 
expect. In the lower panel of Fig. ^ we report the cur- 
rent densities j(x), jc{x), and jr{x) corresponding to the 
electric field and electron-phonon coupling of the mid- 
dle panel. We calculate the current jc by extending the 
single-particle theory of Ref.[l|lto the many-electron case 
within the time dependent Hartree approximation. On 
the time scale of kinetic evolution the current j origi- 
nates from Hamiltonian propagation, while the current 
jc accounts for the charge flow due to inelastic phonon 
scattering and is larger at the contacts where deviations 
from equilibrium are larger. The effect is stronger at the 
injection contact where charge accumulates. The behav- 
ior of jc indicates where local heating should be expected 
to occur. 
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FIG. 2: Two calculated IV characteristics of a DBRTS. The 
curve with the highest peak current (solid line, circles) is the 
one corresponding to the lower value of 70. The energy of the 
dashed line on the left corresponds to twice the energy dif- 
ference between the resonant level and the Fermi level. The 
second dashed vertical line is displaced to the right by twice 
the energy band width. According to a simple model of res- 
onance, transmission should occur in the region between the 
two lines. 



Two I — V characteristics, corresponding to two differ- 
ent values of the electron-phonon coupling parameter 70, 
are reported in Fig. [5J The voltage in the figure is the 
drop across the barrier structure, rather than the total 
voltage drop £ ■ L in the simulation cell. The calculations 
simulate an adiabatic sweep from low to high voltage bi- 
ases and viceversa: at each small increment (decrement) 
of the applied electromotive force we compute a new cur- 
rent and voltage drop using the induced potential ob- 
tained in the previous step as input for the self-consistent 
calculation. Interestingly, the characteristics at lower 70 
exhibits an hysteresis loop which is much less pronounced 
in the characteristics at larger 70. A larger electron- 
phonon coupling also shifts the peak of the resonance to 
lower voltages and reduces the peak to valley ratio. The 
hysteresis loop, which is suppressed by strong electron- 
phonon coupling, originates from charging of the reso- 
nant level, as already found in earlier calculations 23,24 . 

Even though it is not our aim here to reproduce a 
specific experiment, it is interesting to note that by ap- 
propriate choice of the parameters and conditions of the 
calculation, we can reproduce important experimentally 
observed features of a DBRTS like bistability and hys- 
teresis effects^. 

In our approach inelastic collisions provide a source 
of non-linearity other than the tunneling structure. To 
better appreciate this effect we report in Fig. [S] the en- 
ergy distribution of /„„, i.e. the diagonal elements of the 
density matrix, for two values of V, a low voltage below 
peak current and a high voltage close to peak current, 
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FIG. 3: The diagonal elements f nn of the density matrix in 
the steady state. Stronger bias leads to a stronger deviation 
from the unperturbed Fermi-Dirac (FD) distribution. 

respectively. We see that at low voltage the population 
is close to the thermal distribution at temperature T, 
whereas at high voltage the population resembles a ther- 
mal distribution at a temperature higher than T. This is 
a real physical effect: at larger biases the electron distri- 
bution deviates more from equilibrium. We also notice 
(see the inset of Fig. |3J) that the most significant devia- 
tion from thermal equilibrium is a small peak in the tail 
of the distribution, in correspondence with the energy of 
the resonant state. This indicates charging of the reso- 
nant level, an effect which is at the origin of the instrinsic 
bistability of a DBRTS2&. 

The effect of inelastic collisions on transport proper- 
ties can be neglected in the weak bias limit (linear re- 
sponse regime In this regime transport in mesoscopic 
and nanoscopic systems is governed by Landauer's formu- 
las for the conduct ance^i, which can be derived within 
the Kubo formalism^. This is usually done by adopt- 
ing open boundary conditions, but recently a derivation 



based on a ring geometry (in the limit L — > oo) has also 
been reported^. In the linear response regime the con- 
ductance can be calculated by applying a sinusoidal field 
and taking the limit u> — > 0, in absence of coupling to 
the lattice. In our kinetic approach we look instead for 
the steady state solution under a constant applied field. 
This requires coupling to the lattice whenever £ is finite, 
and allows us to consider also strongly nonlinear regimes 
like in the DBRTS case discussed above. 



IV. CONCLUSION 

We have presented a kinetic approach to study quan- 
tum transport at the nanoscale, including the effect of 
inelastic collisions. The main physical approximations 
are a perturbation treatment of electron-phonon scatter- 
ing and a Markov approximation in the dynamics. Use 
of simplified Hamiltonians, as in the example discussed 
above, is by no means necessary. A realistic description of 
the electronic structure within Time Dependent Current- 
Density Functional Theorji22i2ii2£ is possible, and work to 
achieve this goal is at an advanced stage of development. 
More accurate models for the phonons and the electron- 
phonon coupling would be considered in the future. Our 
approach is not limited to steady-state phenomena, but 
could be applied to time-dependent transient behavior as 
well. 
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